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Abstract. The glass phase and its quantum analog are prominent challenges of 
current non-equilibrium statistical mechanics and condensed matter physics. As a 
model system to study the transition from classical to quantum glassy dynamics, we 
propose a setup of laser driven three-level atoms trapped in an optical lattice. Tuning 
the strength of the laser driving to the intermediate level allows one to study the 
transition from a classical Kinetically Constrained Model to the coherent regime. For 
strong driving, Rydberg excitations evolve analogously to defects in the One-Spin 
Facilitated Model, a minimal model known to exhibit glassy dynamics. In our setup, 
the constraints result from the interplay between Rydberg interactions and the laser 
detuning from the Rydberg state. The emerging heterogeneous relaxation timescales 
are tuneable over several orders of magnitudes. In the opposite limit of weak driving of 
the intermediate level, we find an effective cluster model which describes the dynamics 
in a reduced subspace of the allowed number and positions of Rydberg excitations. 
This subspace is uniquely determined by the initial state and is characterized by a fixed 
number of clusters of Rydberg excitations. In addition, we investigate the influence of 
random fields on the classical relaxation. We find that the glassy dynamics can relax 
faster in the presence of weak random fields. 
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1. Introduction 

When disordered materials are quenched to low temperatures, the dynamics can become 
non-ergoclic and the system may fall out of equilibrium [1, 2, 3, 4, 5]. Understanding 
this phenomenon, known as the glass transition, is one of the open questions in 
current statistical mechanics. While a complete theory which satisfactorily explains 
all features of glass forming materials is missing, spin models have been proposed to 
understand particular aspects of the glass problem. For example, static properties 
of glasses are captured by Ising models with random interactions, the so-called spin 
glasses [6, 7, 8]. These models have recently gained renewed interest in the context of 
both quantum adiabatic optimization [9, 10] and, in presence of external random fields, 
of the many-body localization phase transition [11, 12, 13, 14], Spin models which 
study dynamical properties of glasses are so-called Kinetically Constrained Models [15] 
(KCMs). In these KCMs, separation of relaxation timescales as well as space-time 
dynamical heterogeneities [5, 16] result from transitions (spin flips) which are allowed 
or forbidden depending on the state of neighboring spins. Quantum effects in KCMs have 
been recently studied in Refs. [17, 18, 19, 20, 21] in ultracold Rydberg atoms. Thanks 
to their strong and long-range interactions, as well as to their long lifetimes [22, 23, 24], 
Rydberg atoms offer a promising framework to study quantum spin models in- and 
out-of-equilibrium [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]. With the considerable 
recent experimental development in atomic, molecular and optical physics, Rydberg 
atoms can in fact be trapped, laser-excited and detected with single particle resolution 
in state-of-the-art experiments [37, 38, 39, 40, 41, 42, 43, 44], 

In this work we propose a model system based on V-shaped three-level atoms, that 
allows one to study the transition between a classical dissipation-induced glass model 
and a fully coherent quantum system. Atoms are resonantly driven from the electronic 
ground-state | g) to a short-lived low-lying excited state \e) and far-off-resonantly driven 
to a metastable Rydberg state |r) (see Fig. 1). Tuning the driving strength Fl e of 
the | g) |e) transition, allows one to study the crossover from classical to quantum 
dynamics of Rydberg excitations. For large fl e , the relaxation of Rydberg excitations 
resembles the one of defects in a classical KCM, the One-Spin Facilitated Model [15, 18] 
(1-SFM). Here, kinetic constraints are a consequence of the interplay between strong 
interactions of Rydberg-excited atoms [45] and the laser detuning from |r). In contrast, 
for vanishing O e , coherences of the | g) -H- |r) transition dominate the dynamics within 
the lifetime of the Rydberg states [37]. In the case of large laser detuning from the 
Rydberg state the system resembles a transverse Ising model for a particular choice of 
the detuning [46, 47]. We present an effective simple model to describe the evolution in 
this regime which is based on a fixed number of Rydberg clusters with variable number 
of Rydberg excitations. 

In addition to the study of the transition from a glass model to a fully coherent 
integrable model, the proposed framework opens up the possibility to investigate the 
interplay between internal dynamical constraints and external disorder arising from 
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random fields. In onr system, static random fields can be implemented with site- 
dependent random detunings from the Rydberg state. Our results indicate that the 
glassy relaxation of Rydberg excitations toward the steady state can be accelerated by 
the presence of weak random fields. 



Figure 1. The setup we have in mind consists of three-level atoms, each one 
schematized with a ground-state \g), a Rydberg state |r) and an intermediate excited 
state \e). Atoms are coherently laser-driven with Rabi frequencies Sl e (fi r ) between 
\g)-\e) (|g)-|r)) and y e is the decay rate due to spontaneous emission from |e). The anti¬ 
blockade condition V = A r , i.e. when the laser detuning A,, with respect to the \g)-\r) 
atomic transition is equal to the interactions between neighbouring Rydberg-excited 
atoms V (see Subsec. 2.2), ensures that whenever an atom is in |r) ( defect , blue sphere), 
an atom next to it, initially in | g) (facilitated atom, red sphere), will also be favoured to 
be excited to |r). An atom in | g) far from a defect ( non-facilitated atom, black sphere) 
is blocked, because its |g)-|r) transition is far-off-resonant (see Subsec. 4.1). The case 
with additional site-dependent random detunings A' ri is studied in Subsec. 4.3. 

The remainder of the paper is organized as follows: In Sec. 2 the level-scheme and 
the model are introduced; in Sec. 3 we give a brief overview of SFMs. Numerical results 
from quantum trajectory simulations are presented in Sec. 4. In Subsec. 4.2 we study 
the coherent regime and the transition to the rate equation limit and in Subsec. 4.3 we 
analyze the case of finite random detunings from the Rydberg state. The validity of the 
assumptions and approximations of Sec. 4 are discussed in Sec. 5. Finally, in Sec. 6 we 
summarize our conclusions. 
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2. Model 

We consider N atoms individually trapped in a one-dimensional (Id) large spacing 
optical lattice [37, 48, 49, 50, 51] or in a magnetic trap array [52, 53]. Atomic motion due 
to tunnelling between adjacent sites is negligible within experimental timescales. Each 
atom is approximated as a three-level system in the V-configuration, as schematically 
depicted in Fig. 1: the electronic ground-state | g) is laser driven to (i) a Rydberg state 
|r) with a small Rabi frequency h2 r , and (ii) to an excited state |e) with a large Rabi 
frequency Vt e . The intermediate low-lying excited state \e) is rapidly decaying to | g) and 
its short lifetime is given by the inverse of the decay rate 7 e . Two atoms in their Rydberg 
states interact via dipole-dipole V{x) = V/x 3 or van der Waals V{x) = V/x 6 potentials, 
depending on the presence or absence of an applied external held [23], respectively. 

We approximate interactions between Rydberg-excited atoms to be nearest-neighbor, 
i.e. V(l) — V and V(x ^ 2) = 0. Results under this assumption are in qualitative 
good agreement with simulations of the full model including long-range interactions. 
Furthermore, we neglect processes of spontaneous emission from the Rydberg state, 
setting the associated decay rate y r to zero. Restrictions on the validity of these 
approximation are discussed in detail in Sec. 5. 

The Hamiltonian in the rotating frame (with H — 1) is 

# = 5^{y(l2>H + H - c -) + y(l3)H + H - c -) + 

i= 1 ^ 

-v N -1 

- (X, + Ah) | r)(r\i \+V J^||r)(r|i ® |r)(r| m J. (1) 

' i =1 

Here, the sum runs over all N lattice sites, while A r and Ah are the homogeneous and 
site-dependent random laser detunings from |r), respectively. The random detuning 
can, for example, be induced by a spatially inhomogeneous ac-Stark shift acting on 
the Rydberg- or ground-state using either a random light pattern or differential light- 
shifts between atoms in different traps, in close analogy to experiments performed in 
the context of Anderson localization [54, 55]. 

2.1. A single-atom in the V-configuration 

Even though a single three-level systems in the V-configuration has been extensively 
studied in literature, both theoretically [56, 57, 58, 59, 60, 61, 62] and experimentally [63, 

64], let us briefly discuss its properties. For a more complete and detailed discussion, 
we refer the reader to Appendix A. 

The single-atom Hamiltonian of the three-level scheme in Eq. (1), where the 
Rydberg-state is far-off-resonantly laser driven (A r . 0) and no random detuning is 
present (Ah = 0, Vi G [1, iV]), is 

H = y(k)(e| + H.c.) + y(|#)(rj + H.c.) - A r |r)(r 


( 2 ) 
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In the Born-Markov approximation, the master equation describing the dissipative 
evolution of the reduced system density matrix p(t) is in the Lindblad form [65] 

pit) = —i[H, p(t )] + y (2 cp{t)c ] - c f cp(f) - p(t)c ] c ), (3) 

where c = \g)(e\ is the quantum jump operator. Assuming <C kl 2 e /^ e and 
A r f2 e , h2 r , 7 e (the reason behind this choice is explained in Appendix A), we imagine 
to measure photons emitted by our three-level atom. A typical photon counting sequence 
registered by a photodetector would then look like an alternation of time intervals with 
many photons detected ( bright periods, which we label with 0) and no photons detected 
(dark periods, which we label with 1). The transition probability (rate) in going from 
0 to 1 is given by 


T^A,.) 


_ jepppH + 4 a ;)_ 

(ll + 2Ji 2 )[16A; + 4A 2 (7 2 - 2S1 2 ) + SJ‘] ’ 


whereas the transition probability from 1 to 0 is 


( 4 ) 


r^ 0 (A r ) 


i6A? + 4A2( 7 2-2 n 2 e ) + nf 


( 5 ) 


The analytical derivation of these formulas and the physical mechanism behind these 
transitions, known as quantum jumps, are given in Appendix A. 


2.2. Many-body model 

In this Subsection we study the many-body model (A r , V ^ 0) with no random 
detunings. In analogy to the single-atom case [Eq. (3)], the many-body master equation 

N 

p(t) = - i[H , p(t)\ + y J^(2 Ci p(t)c\ - c\dp{t) - p{t)c\ci ), (6) 

i =1 

with Cj = \g)(e\i the quantum jump operator at site i, describes the driven-dissipative 
dynamics of the many-body system. Furthermore, let us assume: (i) Vt r <C Ft 2 e j 7 e , as in 
Subsec. 2.1, for well defined quantum jumps to occur and (ii) A r Ft e ,Fl r ,^ e and (iii) 
V = A r , which we call the anti-blockade condition. The many-body transition rates, 
corresponding to the single-particle rates in Eqs. (4) and (5), can be determined as 
follows. In a Id system, due to the assumption V(x ^ 2) = 0, every process at a certain 
site i can be influenced at most by its two nearest-neighbors, i ± 1. The many-body 
rates can be calculated from Eqs. (4) and (5) with A* = A r — £V the effective detuning 
and £ G [0, 2] the number of neighbors in the Rydberg state. Due to the anti-blockade 
condition, there are three possible configurations: A* = A r if there are no neighbors in 
the Rydberg state (£ = 0), A* = 0 if there is only one neighbor in the Rydberg state 
(£ = 1) and A* = — A r if both neighbors are in the Rydberg state (£ = 2). Since in 
Eqs. (4) and (5), A r appears only in even powers, both off-resonant rates £ = 0, 2 are 
equally suppressed with respect to the resonant rate £ = 1, as highlighted in Tab. 1. 
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A* = 0 

A^ = ±3 7e 

A* = ±10 7e 

r«-‘(A;)/ 7s 

3.00 • 10~ 4 

8.80 • nr 6 

7.54 • 10“ 7 

r 1_>0 (A*)/7 e 

9.00 • 10" 4 

7.14 • 10" 7 

5.64 ■ 10~ 9 


Table 1. Calculation of the creation [r 0_>1 (A*)] and destruction [r 1_>0 (A*)] rates of a 
Rydberg excitation (in units of 7 e ), from Eqs. (4) and (5), for A* = 0 (resonant rates) 
and A* = ± 37 e ,± 107 e (far-off-resonant rates). Other parameters are: fl e = 7 e and 
f i r = 0.037 e . The difference between resonant and far-off-resonant rates is of several 
order of magnitudes and increases with |A*|. 

In absence of conditions (i), (ii) and (iii), the strong transition | g) O |e) might 
significantly increase the probability r 0(1 ^bo)(-j = A r ) of (originally) suppressed processes 
of excitation (deexcitation) to (from) |r) to take place. Since rates differ by several order 
of magnitudes (see Tab. 1), the dynamics shows a large separation of timescales (cfr. 
Sec. 4). 

Allowed (l = 1) and suppressed (£ = 0, 2) processes are analogous to those who 
characterize Spin Facilitated Models, as discussed in Sec. 3. 

3. Spin Facilitated Models 

Spin Facilitated Models (SFMs) are a class of KCMs that exhibit dynamical 
heterogeneities and slowing down of relaxation timescales [15]. Even though they can 
be implemented in a variety of lattice geometries and dimensions, from now on we will 
limit ourselves to Id. Let us now consider the energy function of non-interacting spins 
in an external field, 

N 

E = Kj2 n ^ ( 7 ) 

i =1 

introduced by Fredrickson and Andersen in Refs. [67, 68], where rq = 0,1 is a classical 
spin variable. In our setup, rq = 0 (1) corresponds to the i-th atom being in a bright 
(dark) period. The energy function in Eq. (7) does not feature a finite temperature 
phase transition for a field strength K ^ 0. The equilibrium concentration of defects, 
which we define as spins in the rq = 1 state, reads 

1 N 1 

^ = JV £ ”<(*) = (8) 

i =1 

where fd = l/(fcgT), ks is the Boltzmann’s constant and T is the equilibrium 
temperature of the system. Furthermore, for every K > 0, lim^o d eq = 0, in line 
with the intuition that at low T the concentration of defects should be small. From now 
on, without loss of generality, we set ks — 1 and we fix the temperature scale setting 
K = 1. 

Despite the simple ground-state of Eq. (7), the non-equilibrium dynamics towards 
it can be non-trivial. In general, the evolution of a classical system of spins is described 
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by the following rate equation for the probabilities p( n, t ) of being in the many-body 
spin state specified by the vector n = (n i, ...,rq, ..., nv) at time t: 

<>r ( r ( m ^ n )p( m ^)-r( n ^ m )p( n ^)}- (9) 

rnyn 


Here, r(m —>■ n) is the rate of the transition from m to n 7^ m. The equilibrium 
probabilities 7r(m, t) satisfy the detailed balance T(m — >- n) 7r(m, t ) = T(n —>- m) 7r(n, t ) 
with respect to Eq. (7). Under the assumption of local spin-flips and Glauber 
dynamics [69], the transition rate from (ni,..., n*,..., njv) to (ni,..., 1 — n,,...,njv) at 
site i, is equal to 1/(1 + e@ AE ), with A E = 1 — 2rq the energy difference 

between the final and initial spin configurations. Since in our case A E = ±1, from 
Eq. (8) single-particle rates can be written as 

pn^l-n, = (1 _ de J n . + deq(1 _ n .y (IQ) 


Thus, the asymmetry between the ratio of the rates of creation over destruction of a 
defect 


r 0 ^ 1 d eq 

T 1 ^ 0 1-deq’ 


( 11 ) 


Exes the temperature T at which the system of non-interacting spins thermalizes. (Note 
that T 0 ^ 1 —>■ T 1 ^ 0 for T —)■ 00.) 

Up to now, all transitions between different configurations are allowed. The key 
idea behind SFMs is to introduce restrictions in these transitions. The general rule can 
be stated as: a spin flip is allowed only if there are enough defects in the neighbourhood 
that can facilitate the process. In the Id 1-SFM, these so-called kinetic constraints 
require that a spin can flip provided at least one of its two neighbors is a defect. This 
facilitation condition transforms Eq. (10) into 

r ni->i-ni = ( n ._ i + n . +1 ) |(x _ d eq )ni + d eq (l - rif)}. (12) 


The association between rates in Eq. (12) and the many-body rates calculated in 
Subsec. 2.2 is straightforward: n*_ 1 = n i+ 1 = 0 corresponds to ^ = 0 for which no spin 
flip at site i is allowed, while either {n*-1 = 1, nj + 1 = 0} or {rq_i = 0, n l+ i = 1} 
corresponds to £ = 1 and i is facilitated to flip its spin. The only difference arises in 
the case of {ni-\ = ni + 1 = 1}: while rates in the 1-SFM are doubled with respect to 
the case {?Zj_ 1 = 1, n i+ 1 = 0}, for £ = 2, rates are suppressed as much as in the case of 
£ = 0. This will lead to different steady states which we will discuss in Sec. 4. 


4. Numerical results 

4.1. Semi-classical regime 

In the limit of strong driving to the intermediate state Vt e ~ q e , the dynamics of the 
system is dominated by dissipation. As a method to study the system in the presence 
of both coherent driving and dissipation, we use quantum trajectory simulations 
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(cfr. Appendix A). Fig. 2 depicts a typical trajectory of the Rydberg population r a (t ) 
of several atoms [Eq. (A. 16)]. While for A r = V = 0 [panel (a)], the dynamics of 
Rydberg excitations is unconstrained and all atoms are independent from each other, 
for A r = V 7^ 0 [panel (b)] clear spatial correlations emerge and the dynamics of 
Rydberg excitations resembles the dynamics of defects in the Id 1-SFM [15, 70]. This 
dynamics can be intuitively understood from three different species of atoms, which we 
label, in analogy to the language of SFMs, as: 

• defects, that is atoms in |r); 

• facilitated atoms, that is atoms in either | g) or |e) with one defect as nearest- 
neighbor; 

• non-facilitated atoms, that is atoms in either \g) or |e) with no defect as nearest- 
neighbor. 

The dynamics of defects in the Id 1-SFM is known to be mainly diffusive [70]. Let us 
illustrate the diffusion process of defects in more detail. A typical trajectory is shown 
in Fig. 2(b) and marked by a blue ellipsis. A defect can facilitate the excitation of a 
neighboring atom. In this example, the atom to the right is excited. For a short time, 
both are in the excited state but soon the original atom is de-excited: therefore, the 
defect has effectively moved to the right by one site. Another possible move is highlighted 
in Fig. 2(b) with a red ellipsis. In this second example, the left neighbor is excited but 
later on immediately de-excited. Here, the defect did effectively not move. The diffusion 
process can also be understood quantitatively in some limits. In the regime T • Cl (or, 
equivalently, d eq -C 1), r°^ 1 (0)/r 1 ^°(0) ~ d eq ~ e -/3 and the diffusion constant is 
d e q /2 [15]. In the intermediate regime we studied in this work, the diffusion constant 
does not have such a simple analytical expression, since r°^ 1 (0)/r 1 ^°(0) = 0.33. 

In the following, we study the relaxation of the three different species individually. 
Numerical simulations of Wraj stochastically different trajectories are averaged in order 
to approximate the dynamics of the system described by the master equation (6). 
Results plotting the trajectory averaged Rydberg population Rift) (see Appendix A) 
for each atom i are depicted in Fig. 3. Facilitated atoms (red lines) show a dynamical 
behaviour analogous to the one of resonantly driven single atoms, due to the anti¬ 
blockade condition which acts as a facilitating term in the 1-SFM [cfr. Eq. (12)]. In 
fact, the inset of Fig. 3 shows that, in proximity of t — 0, the profiles of facilitated atoms 
are very well approximated by the resonant (A r = 0) single-atom Rydberg population 
Prrit) (see Appendix A), the associated transition rate being r°^ 1 (0). 

Contrarily, the dynamics of defects (blue lines) and non-facilitated atoms (black 
lines) is slowed down with respect to the one of facilitated atoms. At t — 0, the weak 
transition \g) |r) of a defect is far-off-resonant and r 1 ^°(A r ) is orders of magnitude 

suppressed with respect to r 1 ^°(0) [cfr. Tab. 1]. A defect can in fact relax only after one 
of its neighboring facilitated atoms is excited to the Rydberg state. Similarly, at t — 0, 
the weak transition \g) |r) of a non-facilitated atom is far-off-resonant and r°^ 1 (A r ) 

is orders of magnitude suppressed with respect to r°^ 1 (0). The inset of Fig. 3 shows 
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Figure 2. Comparison of a typical trajectory of Rydberg excitations in the 
unconstrained case A r = V = 0 [panel (a)] and facilitated case A r = V = I0y e [panel 
(b)]. Other parameters are: N = 10, fl e = 7e> W = 0.037 e . The initial condition is 
two Rydberg excitations in * = 1 and i = 6, while all other atoms are in | g). Processes 
of facilitation and diffusion (typical of 1-SFMs) and avoided merging (not present in 
1-SFMs) are highlighted in panel (b). 


that profiles of non-facilitated atoms agree with the far-off resonant (A r ^ 0) single 
atom Rydberg population p rr (t ): non-facilitated atoms, in order to be excited to the 
Rydberg state, have to wait for one of their two neighbors to be Rydberg-excited too. 
Note that for larger N we predict the existence of non-facilitated atoms with profiles of 
Rift) even slower than the ones in Fig. 3. This can be regarded as a form of hierarchical 
relaxation [20] , in that the more distant a non-facilitated atom is with respect to a defect 
at t = 0, the more slowed-down the evolution of Ri(t) will be. 

We now consider the concentration of Rydberg excitations 

1 N 

r ( t ) = ( 13 ) 
i— 1 

The dynamics of r(t) for different N and V, starting with the same initial state is 
depicted in Fig. 4. After an initial transient relaxation, r(t) reaches a steady state 
constant value r = lim t ^. 00 r(t), except in the case of the largest system size N = 10 
and smaller interaction strength V = 3j e , where it increases in time with a small but 
clear trend. We determine the origin of this anomalous slow growth of r(t) in the excess 
number of defects created from non-facilitated atoms within (originally suppressed) 
l = 0 processes, which give a negative temperature contribution in the SFM language 
[compare Eq. (8) with Eq. (11)]. The reason why all reverse i = 0 processes of 
destruction of isolated Rydberg excitations are not relevant, is in the asymmetry between 
off-resonant rates r°^ 1 (A r ) > r 1-5-0 (A r ), which increases with |A r |, as depicted in the 
inset of Fig. 9. More quantitatively, the difference 8r(t) between the excess in the 
concentration of Rydberg excitations at time t and the steady state value r can be 
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Figure 3. Evolution of the trajectory averaged Rydberg population R, (t) for 
TVtraj = 500 in a chain of N = 10 sites. Parameters are: fl e = 7 e , f2 r = 0.037 e , 
A r = V = 37e, A' ri = 0. The initial state is: defects (blue lines) at i = 1 and i = 6; 
facilitated atoms (red lines) at * = 2,5,7,10 and non-facilitated atoms (black lines) at 
i = 3,4, 8,9. Exact master equation evolution of the single atom Rydberg population 
p rr {t) are shown for the initial states p rr (0) = 0 (orange dotted line) and p rr { 0) = 1 
(indigo dotted line). The result r > d eq , suggests that our system lacks of a true 
thermalisation to the ground-state of the 1-SFM (see main text). The shaded thick 
black line is the value of d eq = 0.25 according to Eq. (16). Inset: Resonant (orange 
dotted line) and off-resonant (light brown dotted line) profiles p rr {t) of a single atom 
well approximate the short time dynamics of facilitated and non-facilitated atoms, 
respectively. 


estimated as 

N n{ • r°^(A r .) • t 
N 

where N n f is the number of non-facilitated atoms at t = 0, which we approximate 
constant in time and equal to its value in the initial state. In our case, at the final time 
of the trajectory t max , <5r(t max ) ~ 0.06, which well agrees with numerical results plotted 
in Fig. 4, within statistical errors due to finite IV-traj simulated. Contrarily, the difference 
— 5r'(t ) between the lack in the concentration of Rydberg excitations at time t and the 
steady state value r is 

_ *, (t) = 





where N c \ is the number of isolated defects at t, — 0, which we also assume constant in 
time. Within our parameters, Sr'(t m ax ) ~ 0.003, which is negligible with respect to both 
finite N t ra j errors and 5r(t). From Eqs. (14), and (15) it is clear that, for A r sufficiently 
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Figure 4. Evolution of the concentration of Rydberg excitations r(t) for N = (6,8,10), 
at different interaction strengths, V = (3,10, 20)7 e . Parameters are the same as in 
Fig. 3 and the initial state is, for all plots, two defects as far as possible from each 
other: assuming PBCs, for N = 6 the initial largest distance is 3, for N = 8 is 4, and 
for N = 10 is 5. The length of the dotted line agrees with Sr(t = £ max ) of the main 
text. 


large, l = 0 off-resonant processes are inhibited due to the suppression of r°^ 1 (A r ) and 
r 1 ^°(A r ) and consequently r(t) levels-off, reaching a steady state (see A r > 107 e ). All 
far-off-resonant i = 2 processes in Fig. 4 of merging and splitting of Rydberg excitations 
can be limited by a ’clever’ choice of the initial concentration r(0) and configuration of 
Rydberg excitations. This is the case for, e.g., the initial states of one defect or two 
far defects in Fig. 4, in which the initial number N b of blocked atoms is zero. Blocked 
atoms are atoms with both left and right nearest-neighbour in the Rydberg state: the 
choice of such initial state ensures that, most probably, N b = 0 also at later times and 
we can thus limit the conditions that lead to the occurrence of merging processes. 

Plateaus in r(t) at long times in Fig. 4 show that a steady state is always 
characterized by the general result r > d eq , where 


de q 


7e 2 


2(7 2 + ^r 


(16) 


and not r = d eq , as expected from the theory of 1-SFM. (It can be verified that d eq , in 
the limit <C ff e ,7 e , is well approximated by the single-atom steady state Rydberg 
population p rr (t —> oo) calculated in Appendix B. This confirms that the equilibrium 
ground-state of the 1-SFM is trivial, despite its spatially correlated and slowed down 
transient dynamics.) Our system thus does not thermalize in the long-time limit, but 
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Figure 5. Panel (a): Evolution of the concentration of Rydberg excitations r(t) 
for different initial states r(0). Parameters are the same as in Fig. 3. The initial 
configurations, assuming PBCs, are: r(0) = 0.5 (black line), with defects at i = 
1, 2, 5, 7, 8; r(0) = 0.2 (red line), with defects at i = 1 and i = 5; r( 0) = 0.1 (blue line), 
one defect (due to PBCs, its position is irrelevant). For all initial states, the relaxation 
of r(t) does not converge to d eq . The best agreement with the predicted d eq is obtained 
for the lowest initial concentration of Rydberg excitations, where processes of merging 
(see Fig. 2) are less probable to occur. Panel (b): Evolution of the concentration of 
Rydberg excitations for the initial state r(0) = 0.2 of panel (a) with nearest-neighbour 
(NN), van der-Waals (vdW) and dipolar (dip) interaction potentials. For vdW and 
dip potentials, the range of pairwise interactions have been truncated to third-nearest 
neighbour, which we verified to be a good approximation. 


it converges to a different steady state. Nevertheless, the trend in Fig. 4 suggests that 
for larger N, r decreases approaching d eq . 

Furthermore let us note that, due to finite size effects, r crucially depends on 
the initial configuration of Rydberg excitations, as depicted in Fig. 5(a): while the 
best approximation of the Id 1-SFM result is given, as expected, by the lowest initial 
concentration of defects (blue line), for higher initial concentrations (black line), avoided 
merging processes become relevant and other corrections to r(t), similar to those in 
Eqs. (14) and (15), come into play. Contrarily, in the case in which both l = 2 
merging and splitting processes are present and in the limit of large N , our system 
would thermalize for every initial state in the long-time limit. 

We finally mention that, ideally, one could think of simulating the full dynamics of 
the Id 1-SFM introducing a third laser which weakly couples | g) -H- |r) with a small Rabi 
frequency f! r and detuning 2A r from |r). In this configuration, processes of merging of 
Rydberg excitations would be allowed also for the finite N we simulated numerically. 

4-2. Quantum Regime 

In the limit of vanishing driving of the | g) \e) transition, f \ <C 7 e , f2 r , the dynamics, 
up to times t <C 7 e /f2g, is dominated by | g) |r) coherences. Starting from the fully 
coherent limit Ft e = 0, one can study, increasing hl e , the crossover to the classical rate 
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equation limit of Subsec. 4.1, i.e. VL e ~ y e with Q r <C Ft 2 e /^ e . 

In the case of f2 e = 0, the system Hamiltonian is that of N interacting two-level 
atoms 

N (a 'I N ~ 1 

H t i = jN y^)(r|i+H.c.)-A r |r)(r|i i+V ^||r)(r|i®|r)(r| i+ i|.(17) 

i=1 ^ ' i= 1 

The aim here is to study the non-equilibrium dynamics of the average concentration of 
Rydberg excitations r(t). The Hamiltonian in Eq. (17), assuming V = A r and PBCs, 
can be mapped exactly to the so-called transverse field Ising model [71] . We note that the 
time evolution of the magnetization, which in our formalism can be easily related to r(f), 
following a quantum quench within the ordered (ferromagnetic) phase and across the 
phase transition to the disordered (paramagnetic) phase, has been extensively studied 
in Refs. [46] and [47]. There, the authors provide explicit analytical results in the 
thermodynamic limit N —» oo, corroborated by numerical computations showing that 
finite-size effects are negligible. 

In the following, we show that for V = A r and in the limit <C A r , one 
can approximate the dynamics of excitations of the transverse field Ising model [or 
of Eq. (17)], whose Hilbert space scales exponentially (2^), with the dynamics of 
excitations in a simplified model with only polynomial scaling ( N 2 ). The Hilbert space 
of such a simplified model (Fig. 6) is just a subspace of the transverse field Ising model 
total Hilbert space, and is univocally determined by the initial state of the dynamics. 



Figure 6. Scheme of the effective dynamics governed by the Hamiltonian in Eq. (17) 
with open boundary conditions (OBCs) for t <C A r /Qf, given an initial state with 
energy —V (see main text). The binary values 0 and 1 represent here an atom in | g) 
and |r), respectively. All allowed system configurations are coupled with the single 
atom Rabi frequency f2 r /2 (blue arrows). The black (red) triangle refers to the case 
with N = 3 (N = 4) atoms. In each row, the total number of Rydberg excitations 
N r is constant, independently of N. Increasing the number of atoms from N = 3 to 
N = 4, results in adding a diagonal group of N = 4 new states (highlighted in shaded 
red) and adding a 0 in the rightmost site of all N = 3 states (highlighted in shaded 
black). This scheme can be generalized to an arbitrary N iteratively. 

Let us derive this simplified model from an illustrative example. Considering N = 3 
atoms, H t \ can be represented, in the basis {| ggg), \rgr), \ggr), \grg ), | grr), \rgg ), | rrg), \rrr)}, 
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by the matrix 
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~2 

0 

~ 2 ~ 

-V 

/ 


(18) 


If the initial state has energy — V , the dynamics of the system is limited, for times 
t <C A r /f^, to the subspace \rrr),\ggr),\grg),\grr), \rgg),\rrg) of the full Hilbert 
space. The matrix "H/ =3 , enclosed within the dotted lines in the down-right corner of 
FLf =i in Eq. (18), governs the dynamics of the system in the reduced subspace. 

In our example, the subspace is characterized by the presence of one cluster of 
Rydberg excitations with varying size. Here, a cluster is a group of contiguous Rydberg 
excitations without atoms in | g) in between them. The total number of Rydberg 
excitations N r in the cluster changes during the dynamics (see Fig. 6), but the number of 
clusters does not. At this point we remark that, at times t > A r /f2/ second (and higher) 
order perturbative processes in the small parameter Q r /A r will drive the dynamics 
outside the subspace. This discussion can be extended also for a general N and for 
every initial condition. The number of clusters will be a constant of motion and it 
will always be determined by the initial state. For example, in Fig. 7, we compare the 
evolution of r(t) governed by Fif (green line) and FL^f (red line) for A = 9. We restrict 
the dynamics of the system in the subspace with energy — V starting with one Rydberg 
excitation in the center of the open chain and all the remaining atoms in | g). Numerical 
simulations of the full ("H/ =9 ) and effective ("H/ =9 ) dynamics are in excellent agreement 
for t <C A r /Ft 9 . In this time window, the peaks of oscillations of r(t) have a periodicity 
27r/(fi r /2). 

Remarkably, the dimension of Flff scales as N(N+ 1)/2 contrarily to the exponential 
growth 2 n of FL^i- Such quadratic scaling for large N is understood from the the fact 
that the effective dynamics of our system, in that limit, can be interpreted as hopping 
processes of a single particle in a 2d square lattice, which is the well-known 2d tight- 
binding model. We will comment on this in a few lines. 

It is worthwhile first to study properties of FL^f more in detail. Setting A r = V = 0 
without loss of generality (in the reduced subspace of FL^f, V is just a constant energy off¬ 
set) and considering again the example N = 3, the eigenvalues E^ =3 of "H/ =3 , written 
in ascending order, are 


( a/5 Q r Q r A/5fi r \ 

2~’ 2~ ,0,0, y , ^~ y 


(19) 


In Fig. 8, we plot the spectrum of eigenvalues E\f for different N. Two features emerge: 
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Figure 7. Evolution of r{t) for different fi e , N = 9, and OBCs. Hamiltonian 
parameters are: = 0.037 e , V = A r = 37 e . The initial state is one defect in i = 5, 

all other atoms being in | g). The green curve plots results of simulations evolving the 
full Hamiltonian H^ =9 for = 0. The red curve reproduces with high accuracy the 
exact dynamics evolving the reduced matrix H^ =9 for H e = 0 (see discussion in the 
main text). For all Q e 0, -/V tra j = 500 trajectories have been simulated. 


(■ i ) Independently of N, there is always at least one dark state, that is an eigenstate 
of with Eff = 0. We verified that the total number of degenerate dark states 
increases linearly in N. More specifically, for N ^ 3 [N ^ 2] odd [even], it scales with 
(N + l)/2 [N/ 2], In the case of N — 3, the two (non-normalized) dark states can be 
expressed as 

I D N=3 ) = | rrr) - \grg) and \E N= *) = | rgg) - \grg) + | ggr). (20) 

In particular, \E N ) exists for every N, and it is obtained from the linear combination 
of all states consisting of only one Rydberg excitation with alternating plus and minus 
signs, in analogy to \E N=3 ). Due to the linear scaling in N of the degeneracy of dark 
states, and since bulk properties of the system scale quadratically, we can think of dark 
states also as edge states in the Hilbert space of FUf : for example, \E N=3 ) is localized 
at the lower boundary of the scheme depicted in Fig. 6. 

(ii) The eigenvalue energy spectrum, shown in Fig. 8 for different N, is bounded by 
—2 Q r and 2Q r for large N. This is a consequence of the fact that, as already mentioned, 
our model, in the limit N —y oo, can be mapped exactly to the 2d tight binding model. 
Its dispersion relation looks like 

e(q) = -2J{cos (q 1 s) + cos (qgs)}, (21) 

where e(q) is the energy of the particle at quasi-momentum q = f/ib 1 + g 2 b2 and J is the 
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hopping amplitude between adjacent lattice sites with spacing s. bi )2 = (27r/s 2 )ap 2 are 
the vectors which span the reciprocal lattice generated by a x = s(l, 0 ) and a 2 = s( 0 , 1 ) 

[(1,0) and (0,1) is the standard basis in the Euclidean plane]. Furthermore, due to 
PBCs, qi^s = 27rmi j2 /lV and mi )2 G [-N/2, N/2]: from Eq. (21), e(q) has thus a 
maximum width of 8 J. In our case J = Vt r / 2, and the spectrum width of e(q), 4 fl r , 
coincides with the one observed in Fig. 8 for large N. We stress that this mapping, which 
simplifies a d —dimensional many-body problem with nearest-neighbour interactions to 
a (d+ 1)—dimensional single particle problem, is exact only in the limit N —> oo, where 
boundary effects in our model are absent, i.e. PBCs, assumed in Eq. (21), hold. 


20 


15 


N 


10 


5 


Figure 8. Eigenvalues spectrum for = 0.037 e and A r = 0 at different N. In the 
thermodynamic limit N oo, the spectrum is bounded by ±212,., according to the 
two-dimensional tight binding model dispersion relation [see Eq. (21) and the related 
discussion in the main text]. 

Let us finally comment on the transition from the fully coherent limit to the classical 
limit switching on and increasing 12 e , while keeping fixed all other parameters. Even 
for small f2 e = (0.03,0.1)7 e (blue and orange lines, respectively, in Fig. 7), the system 
dynamics deviates from the coherent case, since Vt e > h 2 r , and oscillations of r{t ) are 
damped by dissipative processes due to the coupling to the short-lived intermediate 
level |e). For t 7 e /fi 2 , r{t) is driven to its steady state value r, which is well 

approximated by p rr (t —> 00 ) [see Appendix B], Upon further increase of h2 e , the rate 
equation limit is recovered, in which r(t) increases slowly with time, oscillations are 
completely suppressed and, as already discussed in Subsec. 4.1, r does not coincide with 
p rr {t —> 00 ). 
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4.3. External disorder 

The external disorder, arising from static spatially-random detunings from |r), A' ri / 0 
in Eq. (1), considerably changes the dynamics of the relaxation. In the quantum 
regime, they are responsible of the many-body localization transition [11, 12, 13, 14]. 

Here, we focus on the effects of random detunings on the relaxation of our system 
in the rate equation limit of the model. We randomly choose the values of A' ri 
within the interval [—A, A] according to a uniform distribution. Moreover, we assume 
0 < A <C A r , in order to avoid accidental balancing of the homogeneous and random 
detunings, A' ri ~ A r . Quantum trajectories are simulated according to the following 
procedure: (i) construct iV r nd different sets of random detunings (A' rl ,..., A ' ri ,..., A' rN ) 
to be added to the homogeneous detuning A r ; (ii) for each random field configuration 
(A' rl + A r ,..., A' ri + A r , ..., A' rN -\- A r ), simulate A tra j stochastically independent quantum 
trajectories; (iii) average the desired observable over the total number of trajectories 
fV tra jfVmd- Note that, for sufficiently large lV md , no extra energy is added to the system, 
because [—A, A] is centered around 0. 

Fig. 9 shows that increasing fV rnd (for a sufficiently large fV tra j), all atoms relax 
faster toward the steady state with respect to the case of zero random detuning. This 
tendency can be understood looking at Tab. 2, which summarizes the effect of a finite 
random detuning on each species. 


Facilitated : A* = 0 

Defect/Non-Fac .: A * = A r 

Blocked-. A * = — A r 

0 + A' ri 

0-IAI.il 

Ar + |Ah| 

Ar-|Ah| 

-Ar + |A'„| 

< 

1 

< 

1 

Speed up 

Speed up 

Slow down 

Speed up 

Speed up 

Slow down 


Table 2. Effect of positive (first inner-columns) and negative (second inner-columns) 
random detunings A' ri £ [—A, A] (with A <C A r ) on facilitated atoms, defects , non- 
facilitated atoms and blocked atoms. The attribute ’Slow down’ (’Speed up’) implies 
that the effective rates r" i_>1_7l ’(A* ± |A' r |) [with m = 0,1] are reduced (increased) 
with respect to zero random detuning rates, (A*). 

As depicted in the inset of Fig 9, for A <C A r , rates with a finite random detuning 
F n i—*i— n i (Ah) of facilitated atoms are larger than the corresponding rates with no 
random detuning r ni_5>1_ni (0) (see the two local maxima of both rates within the beige 
area). As a consequence, at each time, the relaxation of facilitated atoms tends to 
be sped up by any finite random detuning, independently from its sign. Contrarily, 
the relaxation of defects/non-facilitated atoms and blocked atoms is either slowed down 
or sped up, depending on the sign and the magnitude of A' ri (see shaded green and 
blue areas, respectively). As a general rule, blocked atoms and defects/non-facilitated 
atoms are affected in the opposite way by the same random detuning: e.g., for A' ri > 0 
defects/non-facilitated atoms are slowed down, whereas blocked atoms are sped up. 

The result of an averaged (over trajectories and random field realizations) faster 
relaxation to the steady state can be understood more deeply as a consequence of the 
breaking of integrability of the transverse field Ising model the Hamiltonian in Eq. (17) 
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Figure 9. Evolution of the trajectory averaged Rydberg population Rft) for different 
finite random detunings from the Rydberg state, AT 0. The initial configuration and 
parameters are the same as in Fig. 3 and A = 0.57 e (see main text). Different colors 
refer to different numbers of random detuning configurations iV rn( j, for a fixed number 
of trajectories simulated in each random configuration, TVtraj = 100. Inset: Semi-log 
plot of single particle rates T rii_> ' 1_rii (A r ) [with n* = 0,1] from Eqs. (5) and (4) as 
a function of A r . The shaded areas delimit the intervals in which the total detuning 
A* ± |AT| ranges for facilitated atoms (beige), defects/non-facilitated atoms (green) 
and blocked atoms (blue). 


can be mapped to. 

We would like to point out that, in the different limit of large random detunings, 
i.e. A ~ A r , the dynamics is contrarily slowed down, in that facilitated atoms have 
much higher probabilities of turning non-facilitated. 

5. Experimental considerations 

In the following, we discuss and validate the various assumptions made in Secs. 2 and 
4. These are: (i) neglecting the decay from the Rydberg state during the relevant 
experimental timescales, (ii) truncating long-range van der Waals or dipolar interactions 
to nearest-neighbor interactions, (iii) assuming atoms in ground- and Rydberg states to 
be trapped simultaneously in the same optical lattice, and (iv) neglecting heating due 
to photons emitted from the \g) |e) transition. 

(i) The lifetime of Rydberg states is determined by both radiative decay 
(spontaneous emission) to lower-lying electronic energy levels and blackbody radiation 
which induces transitions to higher- and lower-lying neighbouring Rydberg states [23, 
24], While for low-lying electronic levels spontaneous decay to the ground-state 
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dominates, for Rydberg states the interaction with thermally occupied infrared photons, 
which is the main contribution in the overall lifetime, leads to strong mechanical 
effects [72], In general, (alkali-metal) Rydberg states have, compared to low-lying 
excited states, exceedingly long lifetimes oc n 3 (oc n 5 ) for low (high) angular momentum 
states. For n ~ 80, these can be of the order of hundreds of /xs [73]. In our setup, the 
rate of decay from the Rydberg state can be neglected when 


7 r < 




( 22 ) 


4A? 7e 2 + (ill - 4A 2 ) 2 

holds [62], implying that stimulated emission processes from |r) have to be much larger 
than spontaneous ones, which is typically the case for Rydberg states with large enough 
principal quantum numbers n and A r = 0. For A r ^ 0, nonetheless, care has to be 
taken since, contrarily to the case of the general requirement for well-defined quantum 
jumps to occur [62] 


tt 2 r < 


4a; 7 ; + (gt - 4A ;)•- 

il + 4A 1 


(23) 


the condition in Eq. (22) becomes stricter with respect to the resonant case A r = 0. 

(ii) In Fig. 5(b) we check, for a certain initial state and no random detunings, 
the reliability of the approximation of truncating Rydberg interactions to nearest- 
neighbour. While for van der Waals potentials there is no significative quantitative 
disagreement, dipolar interactions result to be detrimental, despite the qualitative 
slow growth of the number of Rydberg excitations is recovered in the long time limit. 
Nevertheless, Rydberg atoms at distances of a few /irn, typical of large spacing optical 
lattices [48, 49, 50, 51, 52, 53], are known to interact with van der Waals potentials [23], 
and our approximation is thus not only qualitatively but also quantitatively justified. 

(iii) It is important, in our setup, that atoms are trapped in the same optical lattice 
during both bright and dark periods. While ground-state atoms can be trapped using 
different state-of-the-art techniques, it has been suggested to optically trap Rydberg 
atoms via ponderomotive potentials from spatially modulated light fields [74] or in dc- 
electric or magnetic fields [75]. In general, electronic ground-states and Rydberg states 
will experience different optical trapping potentials due to their different polarizabilities. 
A promising route to trap ground-state atoms and Rydberg-atoms in the same lattice 
is also to use alkaline-earth atoms as suggested in Ref. [76]. 

(iv) Heating due to photons emitted from the strong transition can be considerably 
reduced adding a finite detuning from |e), A e ^ 0, in order to Doppler-cool atoms with 
the same fluorescence laser. Results of quantum trajectory simulations (not shown) 
indicate that the typical slowed relaxation of atoms is qualitatively unaltered, although 
larger trajectory times with respect to the case A e = 0 are needed in order to observe 
the convergence of r(t ) to its steady state value r. 
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6. Conclusions 

We have shown that the non-equilibrium dynamics of Rydberg excitations in three-level 
atoms can be tuned from the classical to the quantum regime with the laser driving to a 
short-lived intermediate excited state. The anti-blockade condition, where the nearest- 
neighbor interaction among Rydberg levels equals the single atom laser detuning from 
the Rydberg state, implements a form of kinetic constraint. This phenomenon resembles 
the assisted processes of creation and destruction of defects in the One-Spin Facilitated 
Model. The absence of merging processes in our system prevents the thermalization 
to the ground-state of the One-Spin Facilitated Model, despite the concentration of 
Rydberg excitations can be, in the long time limit, stabilized to some steady state 
value for sufficiently large interactions and detunings from the Rydberg state. In the 
intermediate regime of interactions, we observed instead a (tuneable) slow increase of 
the concentration of Rydberg excitations, caused by the appearance of processes which 
are not present in the One-Spin Facilitated Model. In this regime, an additional slow 
relaxation timescale emerges. 

In the coherent limit, the dynamics of the resulting interacting two-level system is 
limited to a subspace of the Hilbert space. This subspace is determined by the choice of 
the initial state, whose number of clusters of Rydberg excitations is conserved. However, 
the size of these clusters changes during the evolution. We derived an effective reduced 
Hamiltonian and verified excellent agreement with simulations of the exact dynamics of 
the full Hamiltonian. 

Finally, we have shown that, in the rate equation limit, the competition of the 
internal disorder of glasses with an externally imposed disorder can result in the slowing 
down but also in the speeding up of the overall relaxation towards the steady state. 

The setup we presented allows for the investigation of the quantum dynamics of Spin 
Facilitated Models in the presence of random fields. As an outlook, we plan to study 
the dynamics of the coherent cluster model in the spirit of many-body localization. 

Acknowledgments 

We thank P. Zollcr for initial discussions which stimulated this work. We further 
acknowledge M. Fleischauer, J. P. Garrahan, I. Lesanovski, B. Olmos, M. Dalmonte, M. 

Heyl and H. Pichlcr. Numerical simulations were performed using QuTip libraries [77]. 

This work is supported by the Austrian Science Fund (FWF): P-25454 N27, the 
SFB FoQuS (FWF Project No. F4016-N23), the European Research Council (ERC) 
Synergy Grant UQUAM, the Austrian Ministry of Science BMWF as part of the 
Unilnfrastrukturprogramm of the Focal Point Scientific Computing at the University 
of Innsbruck. M.M. acknowledges the EU Marie-Curie Program ITN COHERENCE 
FP7-PEOPLE-2010-ITN-265031 for financial support. 


From classical to quantum non-equilibrium dynamics of Rydberg excitations in optical lattices21 

Appendix A. Review of the properties of V-shaped three-level atom 

In this Appendix we review properties of a single atom in the V-configuration and 
describe both quantitatively and qualitatively the physical mechanism of quantum jumps 
in such a setup. 

Let us imagine to measure photons emitted by the three-level atom described by 
the Hamiltonian in Eq. (2) of the main text with a photodetector, which we assume 
to be 100% efficient. The theory of continuos measurement [62] ensures that, under 
continuous monitoring of the strong transition | g) |e) with our photodetector and 
during time intervals where no photons are detected, the system wave function \if(t)) 
evolves according to the Schrodinger equation i-^\if(t)) = H c s\if(t)), with 

H eS = H- ^ct c =H- ^|e)(e| (A.l) 

an effective non-Hermitian Hamiltonian and H given by Eq. (2) in the main text. The 
general solution for the system wave function is 

3 

c n e~ lXnt \u n ), (A.2) 

n= 1 

where X n and \u n ) are the eigenvalues and associated eigenvectors of H c g, 
respectively [66]. The coefficients c n can be determined from the initial condition, e.g. 

1-0(0)) = | g). In the limit <C O e , 7 e , A r , X n and \u n ) can be computed in perturbation 
theory. 


Appendix A.l. Brillouin-Wigner perturbation theory calculation of 
eigenvalues and eigenvectors of H e s 


In the following, we summarize the steps for the calculation of X n and 
\u n ) of H e ff in Eq. (A.l) within a Brillouin-Wigner perturbation theory 
approach. In particular, we focus on the calculation of A e and |« 3 ) (the other 
eigenvalues and eigenvectors can be straightforardly computed following the 
same procedure). In the (|g), |e), |r)} basis, H e g has the following matrix 
form 


H eS = | % 


dg uli'T \ 

~2 ~2 

He 
2 


o 


^ 0 -A r 


Assuming <C Q e ,'y e ,A r , one can split Eq. (A.3) in 

/ 0 ^ 0 \ / 


Hn = 


V 


Oe 

2 

0 


»7e 

2 

0 


0 

-A, 


and Hi = 




\ 


0 

0 

Cl r 

2 


0 

0 

0 


Or \ 
2 

0 
0 


(A.3) 


, (A.4) 
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such that H e ff = Ho + Hi, where Ho is the unperturbed Hamiltonian 
and H\ is the perturbation. Let us introduce a unitary transformation 
{ (|e)(9l-l9>( e l)} w ] 1 i c ] 1 diagonalizes H 0l 


tt t; arctan -j 

U = e 2 1 T'e 


/ Er 


Ho = UHoU^ = 


0 

0 


0 

E n 


0 

0 

-A, 


\ 


7 


(A.5) 


where W is the Hermitian conjugate of U, E m = | y e + y/'y | 4£l\j 

(le + \/7 1 


and E p = 




Hi, Hi = UHiU\ we get 

/ 0 


Applying the same transformation on 


H i = 


Qc \ 
2 


V 


o 

pc 

2 


p) 

. 2 
2 U 


(A.6) 


/ 


where ff c = tt r cosh {|arctanhA^} and Q s = tt r sinh {|arctanh—}. The 
unitary transformation U defines a change of basis from {(#), |e), |r)} to, 
say, {|-), |+), | r)} 




L r ■ 

a 




\+) = i7f\9) + 7f\e) 

u L j- u u y 

l r ) = l r )- 


(A.7) 


We now define two operators V and Q which project onto |r) and its 
complementary subspace, respectively. Since U does not affect the |r) 
subspace [see Eq. (A.7)], V has the same matrix form in either the 
{|0>, |e), |r)} or the {|-), |+), |r)} basis 


/ 0 0 0 \ 


V = 


\ 


0 0 0 
0 0 1 


(A.8) 


while Q can be easily obtained from the orthogonality property of 
projectors, i.e. V + Q = 1, where 1 is the identity operator. The first- 
(E i) and second-order (E 2 ) corrections to the eigenvalue E 0 = [VH 0 V] 3j3 = 
[PH 0 V] 3j3 = —A r . of the transformed unperturbed (H 0 ) Hamiltonian are 


and 


Ei = [PHiV] 3,3 = [VH{P} 3 , 3 = 0 


E 2 = [VHiQH 0 QHiV} 3 , 3 = 


(7e + 2iA r )Q% 


4A r (7 e + 2iA r ) — 2 iFl\ 


(A.9) 


(A-10) 
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Here, [A] t j indicates the i-th row and j-th column of the matrix A and 
IZq = (Eft — QHoQy 1 is the so-called resolvent. Up to second order in 
fl r , A 3 can be approximated with 


A 3 — E n ~ Eq + E\ + E 2 — 

n =0 

_ _ H;j(i 7e - 2A r ) 

r 8 A 2 - 2Q2 -4f 7e A r ' 


(A.ll) 


We can now calculate the first-order correction to |r) originated by the other 
two eigenstates of H 0 , |—) and |+), both mixed to |r) by the transformed 
perturbation Hamiltonian Hi. For example, the eigenvector 1 ^ 3 ) associated 
to A 3 , is: 


l“3> 


(~\Hi\r) (+\Hi\r) 

-A r - E m + —A r - E p " 
Ft r (i' 7e 2 A r ) 

2f 7 e A r - 4A2 + fll ]9) 




|e). 


(A.12) 


2 i 7 e A r — AAj + Ft 2 e 
As mentioned above, the coefficients c n can be calculated from the initial 
condition. For example, setting IV’(O)) = | g), we get 


Q r (i 7e — 2A r ) 
4A^ — fl'l — 2ij e A r 


(A.13) 


The eigenvalues A n are complex numbers all with negative imaginary parts lm{A n } 
[see, e.g., Eq. (A.ll)]. As a consequence, \if(t)) does not conserve the norm during 
the evolution. For -C Hg/ 7e , one of the eigenvalues of H cf[ , A 3 in Eq. (A.ll), 

has a much less negative imaginary part with respect to the other two, |Ina{A 3 }| <C 
|Im{Ai}|, |Im{A 2 }| [ 66 ]. This determines long tails in the delay function 

3 

D 0 (t) = ^|c n | 2 e 2tIm{A " } , (A. 14) 

n= 1 

defined as the conditional probability that no photon has been detected until time t, 
given a photon count has occurred at t, — 0. Note that Do(t) is exactly the norm 
of \if(t)), and its monotonous decay is dominated, in the long time limit, by Im{A 3 }. 
Thus ImjAs} is responsible of the existence of a small but finite probability of a long 
time window with no photon counts. When such a long interval of no photon counts 
occurs, the atom is prepared in the Rydberg state and | f>(t)) is dominated by |« 3 ), 
which in turn can be expressed as |r) plus a perturbatively small admixture of |g) and 
|e) (see Appendix A.l). 








From classical to quantum non-equilibrium dynamics of Rydberg excitations in optical lattices2A 


The wavefunction of trajectories, or system realizations, | ip a {t)) (where a labels 
each of the jV tra j stochastically independent trajectories), can be simulated following 
the procedure below [58]: 

(i) Prepare the system in a certain (normalized) state |'0(^o))a a t the initial time t 0 . 

(ii) Compute the probability of the occurrence of a quantum jump at t 0 from p a (t 0 ) = 
7 e(ifi(to)\c^c\if(to))adt -C 1, where dt is an ideally infinitesimal time step. The 
corresponding probability of no quantum jump is 1 — p a (to) ~ 1. 

(iii) According to these probabilities, randomly choose a no jump/jump event. 

(iv) In case of no jump, propagate \if(t 0 )) a with H eS to obtain the normalised 
wavefunction 1 7p(t 0 + dt)) a = e~ lHcsdt \fj(t 0 )) a /\\...\\ = (1 - iH eS dt)\'ip(t 0 )} a /\\...\\. 
Otherwise, project the wave function into the ground-state, \if(t Q + dt)) a = 

V7ec|^(t 0 ))a/Il-H = I g)- 

(v) Repeat the procedure from (i) to (iv), until the desired final trajectory time f max is 
reached. 


A typical experimental sequence of photon detections looks like an alternation of 
bright and dark periods. During a bright period, labelled by 0, the atomic population 
is rapidly cycled between | g) and \e) and the time between two successive photon 
detections is of the order of the lifetime of |e), q" 1 . In contrast, during a dark period, 
labelled with 1, the atomic population is shelved in |r) and no photons are detected. 
An approximation of the exact time evolution of p(t) governed by Eq. (3) in the main 
text can be obtained simulating A tra j trajectories and finally averaging over them, i.e. 

p(t) = lim ( A - 15 ) 

N t raj^oo A tra j f 

J a=l 

In the following, the observable we are interested in is the Rydberg population, that is, 
the expectation value of the Rydberg projector |r)(r| in | if{t)) a , 


r a {t) = 11 ('r|V^(t)) CK || 2 


(A.16) 


Averaging r a {t) over N tr . A] leads to 

IVtraj 

R(t) = 

JVtra J a=l 


(A.17) 


In analogy to Eq. (A.15), limjv traj ->.oo R(t) = p rr (t), where p rr (t) = (r\p(t)\r) is the 
Rydberg population whose evolution is governed by the exact master equation in Eq. (3). 
Fig. Al shows that R(t) converges to p rr (t) by increasing AT tra j. In the following we 
illustrate three independent methods for the estimation of the average duration of bright 
and dark periods T 0 and Tj, respectively, or equivalently, of the transition probabilities 
(rates) from 0 to 1, T 0 ^ 1 = T 0 -1 , and from 1 to 0, r 1 ^ 0 = Tf 1 . In Ref. [60], these rates 
are found to be 


jpO—>-1 


dp rr 

dt 


1 

t=0 


(A.18) 
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Figure Al. Rydberg population R(t) averaged over A^ tra j trajectories compared to the 
exact evolution p rr {t) governed by Eq. (3). Parameters are: = 7e and fl r = 0.037 e . 

Increasing At ra j, R{t) converges to p rr {t), as expected. 


together with p rr (0) = 0, and 


jpl—>•() 


dp rr (t^j 
dt 


t =0 


(A.19) 


together with p rr ( 0) = 1. Rates can thus be evaluated by linearization of the short-time 
evolution of p rr (t) with the proper initial condition (e.g., from Fig. Al, we can estimate 

r 1 ^ 0 ). 


For comparison, we show the analytical derivation of the rates in Eqs. (4) and (5) of 
the main text. In the long time limit D 0 (t) ~ pe _t//Tl , where p is the (small) probability 
to be in a dark period [61, 62], Comparing it with Eq. (A. 14) one gets 

7 eW 


r i_5>u (A r ) = —2Im{A 3 } = 


16A4 + 4A r 2 ( 7e 2 -2Q2) + n4> 


(A.20) 


and 

p 

where both results are obtained substituting A3 and C3 calculated in Appendix A.l. In 
order to compute r 0 ^ 1 (A r ), we note that, during a bright period, |r) is almost never 
populated and the effective emission rate will be that of the (resonantly driven) two- 
level system | g) |e), i.e. 7 e fl 2 /(7 2 + 2hl 2 ). Then, r°^' 1 (A r ) is simply the two-level 
emission rate multiplied by p [66], 






(7 


e + 4A 2 ) 


16A4 + 4A 2 ( 7 2 -20 2 ) + ^’ 


(A.21) 


r^Ar) 


_7e^(7 e 2 +4A r 2 )_ 

(7 2 + 2W e )[16At + 4A 2 ( 7 2 - 2kll) + ' 


(A.22) 
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Eqs. (A.20) and (A.22) are exactly Eqs. (4) and Eqs. (5) in the main text, respectively. 
If A r = 0 they simplify to 


■' l —>0 j _ j -' l —>-0 


7 


Q 2 

L r 




and 


''O —y 1 ^ q ^ _ j^O—>1 _ 


7 e 3 ^ 




.(A.23) 


The third independent method relies on the direct numerical calculation of the 
average over trajectories of the time-averaged length of dark and bright periods in each 
trajectory. Results for T 1 ^ 0 , obtained with the three different methods, are compared 
in Tab. Al: 



ME 

PT 

QT 

r^°/7e 

9.02 • 10~ 4 

9.00- 10" 4 

9.13 • 10~ 4 


Table Al. Estimates of the resonant [A r = 0] rate r 1 ^ 0 (in units of 7e ), according 
to: the short time slope of p rr (t ), obtained as the exact solution of the master equation 
(ME) [Eq. (3) of the main text], perturbation theory (PT) [Eq. (A.23)] and quantum 
trajectory simulations ( QT ). In all three methods, parameters are: fl e = 7e and 
Pl r = 0.03 7e . Within the QT method, simulations of 7V tra j = 6000 trajectories up to 
times 7 e t = 10® have been performed. 


Appendix B. Steady state Rydberg population 


The steady state solution of the single atom master equation can be obtained setting 
p(t) = 0 in Eq. (3) of the main text. This is equivalent to calculate the kernel of the 
matrix of the corresponding optical Bloch equations 
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P, (B.l) 


where p = (p eg , p ge , Pee, P gg , p rg , p gr , Prr, Per, Pre ) T , p a p = (a\p\/3) and a, (3 = g,e,r. 
(Here, we have suppressed the explicit time dependence of all terms, i.e. p a g = p a p(t), in 
order to improve the readability.) Since the matrix in Eq. (B.l) is singular, we have an 
infinite number of steady state solutions parametrized by the (non-normalized) vector 


f . 7e • 7e 

ysv _! svsv 


7 1 + + n? 




, 0 , 0 , 


y+ 



T 


(B.2) 


pit —» oo) 
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Using the property of trace conservation of p(t), Tr{p(t)} = p gg (t) + p ee (t ) + p rr (t) = 1, 
the normalized steady state Rydberg population can finally be written as 


p rr {t —>• oo) 


7 e 2 + ^ 


2 

r 


2 ( 7 i + ni + n-i)' 


(B.3) 


Notice that Eq. (B.3) converges, for <C fi e ,7e, to the equilibrium concentration of 
Rydberg excitations d eq [Eq. (16) in the main text]. 
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